!
! Generates initial condition for running mhd2d code
! Needs mhd-ic.inp for parameters specification
! Also writes there some diagnostics
!
!=====================================================================
	MODULE size_params
		INTEGER ::   nx, ny
		INTEGER ::   nx2, nx1, nxx, ny2
		INTEGER ::   lx,ly,lxx
		INTEGER ::   n2d,n2dh 
		SAVE
	END MODULE size_params
!=====================================================================
	MODULE fftw_plans
		!For FFTw to work correctly, the size of these integers must
		!match the size of pointers in the language c.
		INTEGER, PARAMETER   ::   ikind = 8       !.. For alpha/64-bit machines
		!INTEGER, PARAMETER   ::   ikind = 4       !.. For most 32 bit machines
		!INTEGER, PARAMETER   ::   ikind = 7       !.. For g77 use (eg, PSR)
		INTEGER( kind=ikind) ::   plan1Df, plan1Db, plan2DRC, plan2DCR
		!INTEGER              ::   plan1Df, plan1Db, plan2DRC, plan2DCR
		SAVE
	END MODULE fftw_plans
!=====================================================================
	PROGRAM ic

	USE size_params
	USE fftw_plans

	IMPLICIT NONE
	INCLUDE 'fftw_include.h'

	COMPLEX, ALLOCATABLE, DIMENSION (:,:)     ::  wk,wh
	REAL, ALLOCATABLE, DIMENSION (:,:)        ::  k2
	REAL, ALLOCATABLE, DIMENSION (:)          ::  kx,ky

	REAL :: time0, efinal, zk0v, qv, zk, zranv1, zranv2, wr, wi
	REAL :: ev, ens, pals, fact, gasdev
	INTEGER :: k_1, k_2, iseed, i, j, iwk
	integer :: j1, j2 !ale

!------------------------------------------------------------
	OPEN(16,file='hd-ic.inp')
	read(16,*) nx
	read(16,*) ny
	read(16,*) k_1
	read(16,*) k_2
	read(16,*) efinal
	read(16,*) zk0v
	read(16,*) qv
	read(16,*) iseed
	read(16,*) iwk

	nx2=nx/2
	nx1=nx2+1
	nxx=nx+2
	ny2=ny/2

	lx=nx2
	ly=ny
	!lxx=lx
	lxx=lx+1

	ALLOCATE( wk(lxx,ny) )
	ALLOCATE( k2(lxx,ny), kx(lxx), ky(ny) )

	IF (iwk.eq.0) THEN
		ALLOCATE( wh(lxx,ny) )

		CALL set_plan2D (nx,ny)                ! FFTW plans
		n2d = ny * (nx + 2)
		n2dh = n2d / 2

	ENDIF

	CALL init_kays (kx, ky, k2)
	time0 = 0.0
	CALL srand(iseed)

	!Generate fourier coefficients wk
	wk = 0.0

	!ale
	j1 = int(ny/4)
	j2 = int(3*ny/4)
	DO i=1, lxx
		zk = sqrt(k2(i,j1))
		if ((zk.lt.float(k_1)).or.(zk.gt.float(k_2))) go to 2
		wk(i, j1) = cmplx(1.0, 0.0)

2		zk = sqrt(k2(i,j2))
		if ((zk.lt.float(k_1)).or.(zk.gt.float(k_2))) go to 1
		wk(i, j2) = cmplx(-1.0, 0.0)

1		CONTINUE
	ENDDO

	!version de dmitruk de la definicion de los wk
! 	DO j=1, ny
! 		DO i=1, lxx
! 			zk = sqrt(k2(i,j))
! 			if ((zk.lt.float(k_1)).or.(zk.gt.float(k_2))) go to 1
! 
! 			zranv1 = gasdev(iseed)
! 			!... zran is gaussian with zero mean and unit variance
! 			wi = zranv1 * zk / sqrt(1. + (zk/zk0v)**qv )
! 			zranv2 = gasdev(iseed)
! 			wr = zranv2 * zk / sqrt(1. + (zk/zk0v)**qv )
! 
! 			wk(i,j) = cmplx(wr,wi)
! 
! 1			CONTINUE      
! 		ENDDO
! 	ENDDO

	!Impose reality conditions
	CALL preal2D (wk)
	CALL zero_kyface (wk)
	!CALL zero_kxface (wk)

	!... calculate the energy and globals
	CALL globals(wk,k2,ev,ens,pals)
        
	!... normalize to efinal
	write(16,*)' '
	write(16,*)' COMPUTED INITIAL CONDITION CHARACTERISTICS '
	write(16,*)' unnormalized ev = ',ev
	fact  = sqrt(efinal/ev)

	wk = wk * fact

	!...  calculate the energy and globals again
	CALL globals(wk,k2,ev,ens,pals)

	write(16,*)
	write(16,*)' AFTER NORMALIZATION'
	write(16,*) 'EKIN         =', ev
	write(16,*) 'ENSTROPHY    =', ens
	write(16,*) 'PALINSTROPHY =',pals
	close(16)
      

	!Writes initial condition
	IF (iwk.eq.0) THEN
		wh=wk
		CALL rfftwnd_f77_complex_to_real(plan2DCR,1,wh,1,n2dh,0,0,0)
		OPEN(1,file='w_0000',form='unformatted')
		CALL write_xarr (wh,time0)
		CLOSE(1)
	ELSE
		OPEN(1,file='wk_0000',form='unformatted')
		WRITE(1) time0,wk
		CLOSE(1)
	ENDIF

	END PROGRAM ic
!=====================================================================
	SUBROUTINE set_plan1D (n)

	!Set the FFTW 1D plans
	!Plans are passed to calling programs through the MODULE fftw_plans

	USE fftw_plans
	INCLUDE 'fftw_include.h'
	INTEGER, Intent( in)     ::  n

	CALL fftw_f77_create_plan( plan1Df, n, FFTW_forward,    &
							FFTW_measure+FFTW_in_place)

	CALL fftw_f77_create_plan( plan1Db, n, FFTW_backward,   &
							FFTW_measure+FFTW_in_place)

	END SUBROUTINE set_plan1D
!=====================================================================
	SUBROUTINE set_plan2D (n1,n2)

	!Set the FFTW 2D plans
	!Plans are passed to calling programs through the MODULE fftw_plans

	USE fftw_plans
	INCLUDE 'fftw_include.h'
	INTEGER, Intent( in)  ::  n1, n2

	CALL rfftw2d_f77_create_plan( plan2DRC, n1, n2,       &
				FFTW_real_to_complex, FFTW_measure+FFTW_in_place)

	CALL rfftw2d_f77_create_plan( plan2DCR, n1, n2,       &
				FFTW_complex_to_real, FFTW_measure+FFTW_in_place)


	END SUBROUTINE set_plan2D
!=====================================================================
	SUBROUTINE write_xarr (aa,time)

	USE size_params
	IMPLICIT NONE
	REAL, Intent(in)           ::   aa(2*nx1,ny)
	REAL, Intent(in)           ::   time
	INTEGER                    ::   i,j

	write(1) time,((aa(i,j),i=1,nx),j=1,ny)

	RETURN
	END SUBROUTINE write_xarr
!=====================================================================
	SUBROUTINE init_kays (kx, ky, k2)

!-------------------------------------------------------
!     kx = 0,1,2,...,nx/2  distributed on proc
!     ky = 0,1,2,...,ny/2,-(ny/2-1),...,-2,-1
!     k2 = kx*kx + ky*ky
!
!     Adjusted so that largest magn ky is POSITIVE
!-------------------------------------------------------

	USE size_params

	REAL, Intent( out), DIMENSION (lxx)        ::  kx
	REAL, Intent( out), DIMENSION (ny)         ::  ky
	REAL, Intent( out), DIMENSION (lxx,ny)     ::  k2
	INTEGER ::  i, j

	do i = 1, lxx
		kx (i) =  float(i-1)
	end do

	do j = 1, ny
		ky (j) =  float( MOD(j-1+ny/2,ny) - ny/2)
	end do
	ky (ny/2+1) =  float(ny/2)                       !.. Make ky_max > 0

	do j = 1, ny
	do i = 1, lxx
		k2 (i,j) = kx(i)*kx(i) + ky(j)*ky(j)
	end do
	end do

	END SUBROUTINE init_kays
!====================================================================
	SUBROUTINE psik (wk,fk,k2)

	USE size_params
	COMPLEX, Intent(in)      ::  wk(lxx,ny)
	COMPLEX, Intent(out)     ::  fk(lxx,ny)
	REAL, Intent(in)         ::  k2(lxx,ny)
	INTEGER                  ::  i, j

	do j=1,ny
		do i=1,lxx
		if (k2(i,j).ne.0.) fk(i,j)=wk(i,j)/k2(i,j)
		end do
	end do
	if (k2(1,1).eq.0.) fk(1,1)=0.0

	RETURN
	END SUBROUTINE psik
!====================================================================
	SUBROUTINE preal2D (aa)
	!for use with (kx,ky) MHD2D code
	!Imposes the interior AND edge reality conditions for a
	!pseudospectral algorithm.

	USE size_params
	IMPLICIT NONE

	COMPLEX,    Intent( inout)  :: aa(lxx,ny)
	INTEGER                     :: jy, jcjg

	!Modes (0, ky) are conj to (0, -ky)
	DO  jy = 2, ny/2
		jcjg = ny - jy + 2
		aa (1,jcjg) = CONJG ( aa (1,jy) )
	ENDDO

	!Do the special point k=(0,0) explicitly [=0, not just pure real]:
	aa(1,1) =  0.0

	RETURN
	END subroutine preal2D
!=====================================================================
	SUBROUTINE zero_kyface (aa)

	USE size_params
	IMPLICIT NONE

	COMPLEX,    Intent( inout)  :: aa(lxx,ny)

	aa (:, ny/2+1)  =  0.0

	RETURN
	END subroutine zero_kyface
!=====================================================================
	SUBROUTINE zero_kxface (aa)

	USE size_params
	IMPLICIT NONE

	COMPLEX,    Intent( inout)  :: aa(lxx,ny)

	aa (lxx, :)  =  0.0

	RETURN
	END subroutine zero_kxface
!=====================================================================
	SUBROUTINE globals(wk,k2,evt,enst,palst)

	USE size_params
	USE hd_params
	IMPLICIT NONE

	COMPLEX, Intent(in)      ::  wk(lxx,ny)
	REAL, Intent(in)         ::  k2(lxx,ny)
	REAL, Intent(out)        ::  evt,enst,palst
	INTEGER                  ::  i,j
	REAL :: wr, wi, ks

	!evt   = kinetic energy = 1/2 sum (wk^2 / k^2)
	!enst  = enstrophy = 1/2 sum (wk^2)
	!palst = palinstrophy = 1/2 sum (k^2 * wk^2)

	evt   =  0.0
	enst  =  0.0
	palst =  0.0

	DO j=1,ny
		DO i=1,lxx
			!Count only the energy of the mode if ( kx > 0 .or. ky > 0 )
			if ( i.gt.1 .or. (j.gt.1 .and. j.le.ny2) ) then !ale: no entiendo esta linea.
				wr = real(wk(i,j))
				wi = aimag(wk(i,j))
				ks = k2(i,j)

				evt   =  evt   +  (wr*wr + wi*wi)/ks
				enst  =  enst  +  (wr*wr + wi*wi)
				palst =  palst +  (wr*wr + wi*wi)*ks
			ENDIF
		END DO
	END DO

	RETURN
	END SUBROUTINE globals
!=====================================================================
	FUNCTION GASDEV(IDUM)
	DATA ISET/0/
	common /gsave/gset

	IF (ISET.EQ.0) THEN
1		V1=2.*RAND()-1.
		V2=2.*RAND()-1.
		R=V1**2+V2**2
		IF(R.GE.1.)GO TO 1
		FAC=SQRT(-2.*LOG(R)/R)
		GSET=V1*FAC
		GASDEV=V2*FAC
		ISET=1
	ELSE
		GASDEV=GSET
		ISET=0
	ENDIF

	RETURN
	END
!=====================================================================
